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An accepted model reduction technique is balanced 
cation, by which negligible states of a linear system of 
are determined by balancing the systems controllability and 
observability gramian matrices. To be applicable for nonlin- 
ear system this method was enhanced through the empirical 
gramians, while the cross gramian conjoined both gramians into 
one gramian matrix. This work introduces the empirical cross 
gramian for square Multiple-Input-Multiple-Output systems as 
well as the (empirical) joint gramian. Based on the cross gramian, 
the joint gramian determines, in addition to the Hankel singular 
values, the parameter identifiability allowing a combined model 
reduction, concurrently reducing state and parameter spaces. 
Furthermore, a controllability and an observability based com- 
bined reduction method are presented and the usage of empir- 
ical gramians is extended to parameter reduction in (Bayesian) 
inverse problems. All methods presented are evaluated by nu- 
merical experiments. 
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1 Introduction 



Efficient reduction of large-scale nonlinear control systems is a challenging task. 
Even more in the case of parametrized systems, in an inverse problem set- 
ting, where a combined reduction of parameters and states is targeted. For 
instance, a network with many nodes and unknown connectivity, modeled as 
a parametrized nonlinear dynamic system, requires much time during param- 
eter estimation due to system size and parameter count. With the aim to 
lower computational complexity, the parameter and state spaces are confined 
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to low-dimensional subspaces. One class of model order reduction techniques is 
concerned with determining projections to such subspaces, mapping the high- 
dimensional model to a reduced surrogate model. Here, the projection based 
method of balanced truncation is the origin of the proposed combined reduction 
of state and parameter space. 

Since the number of a systems inputs and outputs usually remains fixed, the 
maps to and from the intermediary states characterize the reducibility of a 
system ([T]). The balanced truncation approach, introduced in [2], categorizes 
the states of a system in terms of controllability and observability. From the 
input-to-state map the systems controllability is computed, quantifying how 
well a state is driven by input. The state-to-output map is used to compute 
the systems observability, quantifying how well changes in a state are reflected 
by the output. Excluding the least controllable and observable states from a 
model motivates a projection from the full to a truncated model. Given a linear 
time-invariant control system. 



x{t) = Ax{t)+Buit) 
y{t) = Cx{t), 




> 



this systems controllability and observability can be assessed through the asso- 
ciated gramian matrices Wc and Wq, which are classically computed as the so- 
lutions of the Lyapunov equations AWc + WcA* = -BB* and AWq + WqA* = 
—C*C. To make a compound statement about controllability and observability, 
Wc and Wq have to be balanced. The singular values of the resulting bal- 
anced gramian correspond to the Hankel singular values of the system, with 
their magnitude describing how controllable and observable the associated state 
is. A description of the balancing procedure can be found, for example, in [5] 
or [1]. 

This work focuses on cross gramian based methods for model reduction. The 
cross gramian was introduced in [5]; it combines controllability and observabil- 
ity information into a single gramian matrix and is elaborately described in [B]. 
The empirical gramians were introduced in the works |7], [S] and enable the 
computation of gramian matrices for nonlinear systems by mere basic matrix 
operations. This concept was extended among others, in [5] allowing nonlinear- 
ity measurement and in [lOj providing more general input signals. Particularly 
noted should be [TT] and [T^] for developing the empirical cross gramian for 
SISO systems in the context of sensitivity analysis. For the gramian based pa- 
rameter reduction, the groundwork has been laid by fl3J , from the observability 
standpoint and by [T3], from the controllability point of view. 
To begin, the cross gramian and its properties are reviewed in section 2) Then 



the empirical cross gramian (for MIMO systems) is developed in section 3 Sec- 
|tion 4| introduces combined state and parameter reduction in two variants. First, 
a controllability and second, an observability based approach; the latter is en- 
hanced to a cross gramian based combined reduction, that is presented in |sec- 



tion 5 A procedure to reduce parameters spaces in Bayesian inverse problems, 
prior to parameter estimation, is suggested in [section 6] Finally, numerical ex- 
periments are conducted in [section 7| comparing the here presented methods for 
linear and nonlinear systems. 
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2 Cross Gramian 



The cross gramian Wx (also known as: Wco) was introduced in a sequence of 
works (0, dg, [Ml, nil, HSI, jH], EO]) and encodes controllability and ob- 
servability into a single matrix, yet it can only be computed for square systems, 
meaning a system with the same number of inputs and outputs. In this case 
the cross gramian is the smallest semi-positive solution of the Sylvester equa- 
tion AWx + WxA — —BC, of which approximate solutions were discussed in 
|21| . [55] and [53] • The cross gramian can also be specified by the product of 
input-to-state and state-to-output maps: 



Wx= I e^^BCe'^^dT, ( (1) 



in case the underlying system is asymptotically stable. If the system is further 
on symmetric, the following relation between the cross, the controllability and 
observability gramian holds ([S]): 

W 

Wl = WcWo \\iWx)\ = VHWcWo). 

And as shown in [55', Lemma 2.1], the eigenvalues of the cross gramian corre- 
spond to the Hankel singular values. 

While a for Single-Input-Single-Output (SISO) system this is always true (L5j), 
a linear Multiple-Input-Multiple-Output (MIMO) system not only requires the 
same number of inputs and outputs, but also the system gain G = —CA~^B 
to be symmetric ([5T]), which is equivalent to the existence of a symmetric 
transformation J, with AJ = jA^ and B — C^J. Trivially for J = 1 the 
system would be restricted hy A — A^ and B — C'^ . 

As presented in [5_ the trace of the cross gramian equals half the gain of a SISO 
system 




tT{Wx) = ~lcA-'B. 



The trace being the sum of eigenvalues, this relationship associates the eigen- 
values with the system gain. This result was used in [11] and [12] for parameter 
identification purposes, using the system gain as a sensitivity measure. For 
MIMO systems, here an extension to ^ Theorem 3] is developed: 

/ 

Corollary 

Given a linear, square MIMO system, the trace of the cross gramian relates to 
the system gain as follows: 



^A^x)^-\'^ACA-^B). 



Proof. 

For an asymptotically stable system, the trace of the cross gramian (in the form 



of equation 1 1 is given by: 

POO 

tr{Wx)=tT{ e^^BCe^^dr) 
Jo 
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tr{e^^BCe'^^)dT, 



and since the trace of a product of compatible matrices is commutative, 

/■OO 

Jo 

POO 

= tr( / Ce^'^^Bdr) 
Jo 

= tr(C / e^^^drB). 





tr{Ci--A-')B) 
-^-tT(CA-^B). 



Employing the cross instead of controllability and observability gramian, means 
only a single gramian has to be computed. And since no balancing is required, 
the truncation procedure can be simplified to direct truncation ([H], |25l Ch. 
7], [HI Ch. 12.3]). The projection is given by the singular value decomposition 
of the cross gramian, 

jr 

UDV = svd{Wx). 

Then U is truncated based upon the singular values D, yielding a projection 
matrix U. Applying U to A,B,C results in the reduced system: 

jir{t) =U^AUxr{t) + U^Bu(t) 
^ ^ly(i) = ClJx,{t). 

Lastly, the cross gramian has the benefit of conveying more information then 
the^controllabihty and observability gramian (jS]). 

3 Empirical Cross Gramian 

For nonlinear models, 

X = f{x,u,9) 
y = 9{x, 0) 
x(0) = xo, 

the above procedure is not viable. In [8J and |26j the concept of empirical 
(controllability and observability) gramians was introduced. This is an extension 
of the PC A method which uses the SVD of a snapshot matrix ([37]). Subse- 
quently this approach and its field of application was advanced by j28], [29], [30] 
and [3T]. The empirical gramians not only correspond to the classic gramians 
for linear systems ([7]), but might even carry more detailed information on the 
underlying system, as shown in |32| . Empirical gramians are based on averaging 
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the output of a system that is perturbed in inputs and initial states. Initially 
only for delta impulse type input, this ansatz was extended in |10| for more 
general input configurations under the name of controllability and observability 
covariance matrices. 

Next, the necessary perturbation sets are defined; these should reflect the oper- 
ating range of the underlying system. E and F are sets of standard directions 
for the inputs and initial states. Sets S and T are orthogonal transformations 
(rotations) to these standard directions of inputs and initial states respectively, 
while Q and R hold scales to these directions: 

E = {e^e R™; ||e,|| = 1; e,e, =0:i = l,. 
F = {f,eR'';\\M^l;Mj=0;i = l,. 
S = {S, G R™^"'; S*S, = 1; ^ = 1, . . . , s] 
T = {T, e R""""; 7;*r, = 1; i = 1, . . . 
Q = {(7i e R; > 0; i = 1, . . . , g}, 
i? = {r^ e R;r, > 0;i = 1,.. .,r}. j 

The empirical cross gramian for SISO systems was introduced in here, as 
a new contribution, this concept will be generalized to square MIMO systems. 
Hence the scope of the cross gramian is extended to nonlinear systems. For 
square MIMO systems the empirical cross gramian is then defined by: 

Definition (Empirical Cross Gramian) ^ ^ 

For sets E, F, S, T. Q. R and the steady state x,u{x)^y{x) the empirical 
cross gramian relating to the output the states x^^^ and the input u{t) — 
CmSiifj6{t) + u as well as the initial state a;(0) = dj-TiCi + x is given by: 

f h=l k=l 1=1 m=l " 

A.T™''^ (<) = (x"''^'(t) - x"''^ ) 
A/'*(i) = (/'^(t)-y'='*). 
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Along the lines of f7^ and [TT] , the equality of the cross gramian and the empirical 
cross gramian of a linear MIMO system is shown next: 



Lemma (Empirical Cross Gramian) 

For any nonempty sets S, T, Q, R the empirical cross gramian Wx of an 
asymptotically stable linear control system is equal to the classic cross gramian. 

Proof. 

For an asymptotically stable linear control system, Ax{t) = x{t) = e^^Bu{t) 



and similarly Ay{t) = y{t) = Ce xq, thus 



= efTfie^^Bc.^S,f,)ffS^iCe^'d,T,e,) 
= dkC^ejTfe^'BCe^'Tie^ 



1 _ _ 



k=l 

= Wx. 




Following ini, a discrete representation of the empirical cross gramian is given 
here too, due to the situation that data is collected at discrete times. 



Definition (Discrete Empirical Cross Gramian) 
For sets E,F,S,T,Q,R the discrete empirical cross gramian relating to the 
output the states a;''^' and the input u{t) = ci,Sjei5{t) + u as well as the 
initial state x{0) = dkTiCi + a; is given by: 

The computational complexity depends largely on the number of scales and 
rotations of perturbations as well as the cost of integration. 

4 Combined State and Parameter Reduction 

In this section two methods for combined model reduction, allowing simulta- 
neous reduction of state and parameter spaces are introduced; an observability 
based and a controllability based ansatz. 

In |13| the identifiability gramian was introduced, which extends the concept 
of observability from states to parameters. Augmenting the states of a given 
system by including its parameters as constant states, the parameters are treated 
hke states: 



y 



x 



f{x,u,e) 





a; 



y = gix,o) 
(0) = ^^(0) 



These augmented states are initialized by their parameter value. The identifia- 
bility of the parameters is obtained through the observability of these augmented 
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parameter-states. From the observability gramian of this augmented system, the 
identifiability gramian Wi can be extracted via the Schur-complement , 



{Wo 






wj 



Wi 

=^Wi = Wp- W:Wo'^Wa. V 

The identifiability is then given as the observability of the parameters. A param- 
eter reducing projection is computed, similarly to the cross gramian approach, 
through a singular value decomposition of the identifiability gramian. Observ- 
ability of the states is encoded in this combined observability gramian, too; it 
can be extracted as the upper- left N x N matrix from Wj. Thus after a iden- 
tifiability gramian based parameter reduction, and computation of the systems 
controllability gramian, a state reducing balanced truncation can be applied to 
achieve combined reduction. 

In [T^ the per parameter controllability gramians were introduced for sensitivity 
analysis purposes, in which the parameters are treated as additional inputs. 
Controllability gramians are computed for the system without any parameter 
inputs (Wc,o)j and for each parameter as sole input to system excluding other 
parameters as well as regular inputs (Wca)- Then a controllability gramian of 
this system is given by the sum of all sub controllability gramians: 



1=0 <^ 



The sensitivity related controllability information (the traces of the associated 
sub-gramian Wc,i>o of each parameter) can then be accumulated in a diagonal 
matrix, which is called the sensitivity gramian (see |3?|): 






tT{Wc,p) 

The singular value decomposition of Ws induces a parameter projection. In 
combination with an observability gramian, a balanced truncation of states sub- 
sequent to a sensitivity gramian based parameter reduction results in a combined 
reduction. 

> ^ 

5 Joint Gramian 



Now, aggregating the computation of controllability and observability not only 
for states as proposed above, but also for identifiability of parameters leads to a 
reduction method requiring a new single gramian. Additionally to extending the 
states with as many (constant) states as there are parameters; inspired by [34J, 
the inputs axe extended to include controllability of the parameters. To be able 
to apply the cross gramian to this extra augmented system, the outputs have 
to be augmented as well, to keep the system square. The augmented inputs 
act untransformed upon the parameters; consequently, to ensure the systems 
symmetry, the output transformation of the augmented outputs is represented 
by an identity function. This leads to the following new gramian matrix: 



Definition (Joint Gramian) 

An extra augmented system, with v = 0, is given by 



f{x,u,0) 



g{x,e) 
t 

x(0) 




y 

m 



The cross gramian of such extra augmented system is called the joint gramian. 

Applying the empirical cross gramian to this setting allows to compute the 
empirical joint gramian. which is composed of: 



(Wx 










Wj = 



Wx being the cross gramian and Wp the parameter-related gramian. The iden- 
tifiability information on the parameters can be extracted by the Schur comple- 
ment: 



W; 



Wp-W:Wx-^Wa, 



with Wj € R^^^ the cross identifiability gramian, which is usually not equal 
to the identifiability gramian presented in [13_, since here not only observability, 
but also controllability is part of the identifiability. As for Wi and Ws the as- 
sociated parameter projection is computed by the singular value decomposition 
of Wj. Since the (empirical) gramians are symmetric, a variant of the Cholesky 
decomposition, based on [35J, can be used to compute the Schur-complement: 

Wj = Wj + \\Wj\\il 

L*L = cho\{Wj 




u , Wx W* 

Wp 



La Li 

L, 



Wj = 

Wj^Wj~\\Wj\\il. 

For a combined reduction of states and parameter no further gramian has to be 
computed and no balancing is required, since the cross gramian Wx for state 
reduction is given by the upper-left N x N matrix of the joint gramian Wj. 



6 Bayesian Inversion 

The empirical gramian method is especially well suited for parameter identifi- 
ability, for example see [36.. Parameter identifiability or reduction for inverse 
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problems is a complicated setting, since the parameter reduction should occur 
before their estimation. Especially for Bayesian inverse problems, empirical 
gramians can be utilized with ease due to prior distribution on the unknown 
parameters, which is presented next: Assuming a Gaussian setting, the pro- 
vided prior knowledge, given by mean and covariance can be employed to define 
the configuration of the parameters steady states and perturbations. Using the 
prior mean, the steady states of the augmented parameter states are set. 



X ^ f{x,u,E[p{9)]) 
y ^ gix,E[pi9)]), 



while the prior covariance defines the perturbation scales, given a sequence 
a^ G [0,1]: 



uence 



{qi = Gi ^jYai[p{dj)]; i = l...m}, 



Re, = {ri = aiJNai[p{9j)]]i = 1 . . . rt}. 



For the actual parameter reduction, after the identifiability has been determined, 
two reduction strategies are available. First, following the truncation procedure, 
by projecting the priors to a low-dimensional subspace and estimating only a 
low number of parameters. Second, by identifiability analysis, using the eigen- 
vectors associated to the eigenvalues of the gramian, the reducible parameters 
are excluded from the estimation procedure and kept at their prior mean value. 
As described above not only parameter reduction, but also combined reduction 
of states and parameters can be applied in this pre-estimation Bayesian inverse 
problem setting. 

^ / 

7 Implementation and Numerical Experiments 

For an efficient implementation the structure of the gramian computation is 
exploited. First, the empirical gramians allow extensive parallelization. Each 
combination of direction, transformation and scale can be processed separately 
yielding a sub-gramian. Second, the assembly of each sub-gramian can be 
comprehensively vectorized, since it consists of vector additions, inner- and 
outer-products. The output gramian is the normalized accumulation of all sub- 
gramians. For further details about the implementation refer to [33) . 
The empirical gramian framework can be found at |http : / /grEunian . de/emgr . m 



and is compatible with Octave and Matlab while the source code, used for the 
following experiments, can be found at http://gramian.de/himpel3a.tgz. 
The empirical cross gramian and empirical joint gramian are verified with a 
(parametrized) linear system, 

yr x{t) = Aex{t) + Bu{t) 

y{t) = Cx{t), 

and tested with a (parametrized) nonlinear system, with the nonlinearity being 
a component- wise applied sigmoid function, here a hyperbolic tangent (see j37J); 

x{t) = Ae tanh(a;(i)) + Bu{t) 
y{t) = Cx{t). 
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In both cases the matrix Ag is fully parametrized, meaning each component of 
Ag is a separate parameter, hence 9 ^ — . Here the linear and nonlinear 
system are chosen to have 1 = — 4 inputs and outputs as well as = 32 
states and thus P = N"^ = 1024 parameters. Regardless of the Hankel singular 
values, state reductions reduce to Nn = 4 states and parameter reductions to 
Pfi = 16 parameters. The parameters constituting Ag are generated randomly, 
but stability and symmetry of the system are ensured. For parameter and 
combined reductions the experiments are treated as Bayesian inversion, reducing 
the parameter space before estimation that is initialized with the prior mean. 
Input matrix B is also random and the output matrix is given by C = to 

fij^ [fijni 1%:^ 



guarantee symmetry. Following, each figure (fig. 1 fig. 2 



fig. 6 1 shows the relative error of the reduced to the full model for each state 
over time for an exemplary system. ^ ^ 

The state reduction using the empirical cross gramian {Wx) is compared to 
balanced truncation of empirical controllability gramian (Wc) and empirical 
observability gramian (Wq), and also to a variant of balanced proper orthogonal 
decomposition from [38^ and [39]. ^ 





offline 


online 


L2-f 




time [s] 


time [s] 


error 


original system 





0.0061 





balanced POD 


0.1279 


0.0112 


0.0139 


balanced truncation {Wc,Wo) 


0.9747 


0.0110 


0.0223 


direct truncation (Wx) 


0.6021 


0.0111 


0.0149 



Table 1: Linear State Reduction, averaged over 10 systems. 





offline 


online 


L2- 




time [s] 


time [s] 


error 


original system 





0.0063 





balanced POD 


0.1330 


0.0125 


0.0154 


balanced truncation {Wc.Wq) 


1.0462 


0.0117 


0.0370 


direct truncation {Wx) 


0.6538 


0.0117 


0.0151 



Table 2: Nonlinear State Reduction, averaged over 10 systems. 



For either, the linear (table 1) and nonlinear (table 2 1 state reduction the cross 
gramian based direct truncation requires less time than the balanced truncation 
using controllability and observability gramian. Naturally, balanced POD is 
faster, but produces a lower error only for the linear case; for the nonlinear 
state reduction the empirical cross gramian reduces with the least error. 
For parameter and combined reduction the parameters constituting Ag are 
treated as unknown and are initialized with flat priors. Thus each online phase 
solves an inverse problem by estimating the parameters using unconstrained 
optimization. The parameter reduction performance, of the empirical joint 
gramian (Wj) is assessed in comparison to the empirical sensitivity gramian 
{Ws, see [33]), and the empirical identifiability gramian {Wj, |13|). 
In all cases (see table 3 and table 4 1 , the online time is reduced by approximately 
two orders of magnitude. Linear and nonlinear parameter reduction is completed 
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Figure 1: Exemplary absolute error comparison of state reduction for a linear 

system by balanced POD, balanced truncation (using empirical controllability 
gramian and empirical observability gramian) and direct truncation (using the 
empirical cross gramian). 
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w 

offline 


online 


L2- 




time [s] 


time [s] 


error 


original system 





353.6360 


0.0228 


t direct truncation (Ws) 


38.2590 


9.3704 


1.5826 


direct truncation (Wj) 


133.5540 


5.9549 


0.3361 


direct truncation (Wj) 


1143.8626 


6.0806 


0.3204 



Table 3: Linear Parameter Reduction, averaged over 10 systems. 





offline 


online 


L2- 




time [s] 


time [s] 


error 


original system 





355.4448 


0.0091 


direct truncation {Ws) 


40.4850 


8.6330 


0.9922 


direct truncation {Wi) 


136.2866 


7.5322 


0.0988 


direct truncation {Wj) 


1163.3324 


7.5328 


0.0578 



Table 4: Nonlinear Parameter Reduction, averaged over 10 systems. 



with the least error by the joint gramian, yet requires significantly more offline 
time than the sensitivity and identifiability gramian. The sensitivity gramian 
Ws demands the least offline time, but the resulting error is up to an order of 
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Figure 2: Exemplary absolute error comparison of state reduction for a non- 
linear system by balanced POD, balanced truncation (using empirical controlla- 
bility gramian and empirical observability gramian) and direct truncation (using 
the empirical cross gramian). 

/ T XT 

magnitude larger compared to the identifiability gramian W7, with its offline 
time still being below the original models online time. 

Finally, the combined state and parameter reduction of the empirical joint 
gramian is benchmarked against the two variants of combined reduction by bal- 
anced truncation described [above] The first variant reduces the parameter space 
by the empirical sensitivity gramian; since a controllability gramian is already 
computed in the process, then with an empirical observability gramian the state 
space is reduced by balanced truncation. The second variant reduces the param- 
eter space by the empirical identifiability gramian, then a state space reduction, 
utilizing balanced truncation of the empirical controllability gramian and ob- 
servability gramian, latter being a byproduct of the identifiability gramian. 





offline 
time [s] 


online 
time [s] 


L2- 
error 


original system 

balanced truncation {Wsi^o) 




39.1512 


353.6360 
12.9656 


0.0228 
0.0611 


balanced truncation (VFc,W7) 


133.6889 


6.0815 


0.1117 


direct truncation (Wj) 


1143.8660 


6.6106 


0.0498 



Table 5: Linear Combined Reduction, averaged over 10 systems. 
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Figure 3: Exemplary absolute error comparison of parameter reduction for a 
linear system by empirical sensitivity gramian, empirical identifiability gramian 
and empirical joint gramian. m 
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^ offline 


online 


L2- 




time [s] 


time [s] 


error 


original system 





355.4448 


0.0091 


balanced truncation {Ws,Wo) 


41.4663 


12.1245 


0.0423 


balanced truncation (VFc, W/) 


136.3902 


8.4331 


0.0438 


direct truncation {Wj) 


1163.3232 


8.0231 


0.0360 



Table 6: Nonlinear Combined Reduction, averaged over 10 systems. 



[Table 5| and [table 6| show, along the the line of state and parameter reduction, 
the joint gramian fitting the original model closest. For either, the linear and 
nonlinear combined reduction, the joint gramian demands the most offline time 
while the online time is comparable to the balanced truncation methods. 

8 Conclusion 

In the case of nonlinear systems, the empirical cross gramian is advantageous 
compared to other gramian based model reduction methods. The cross gramian 
approximates, on average, in state, parameter and combined reduction the sur- 
rogate models with the least error. Moreover, the cross gramian can be ap- 
plied to non-symmetric square systems with broader error bounds ([22l Remark 
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Figure 4: Exemplary absolute error comparison of parameter reduction for 
a nonlinear system by empirical sensitivity gramian, empirical identifiability 
gramian and empirical joint gramian. m 



2.1]). Surrogate models based on the joint gramian (parameter or combined 
reduction), though the computationally most costly, provide the best results in 
combined (state and parameter) reduction. Due to the longer assembly times of 
the joint gramian it is suggested to use it only if the dimension of the parameter 
space is of the same scale as the dimension of the state space. Further research 
has to be conducted on applying the empirical cross gramian to nonlinear or 
non-symmetric systems. Either by heuristically applying it f|23j). extending the 
symmetry constraint to nonlinear systems (|40]), or symmetrization (|21|). 
Combined reduction of states and parameters using empirical gramians with 
either sensitivity and observability gramian, controllability and identifiability 
gramian or the joint gramian is applicable to (Bayesian) inverse problems. While 
the combination of sensitivity and observability gramian has the fastest offline 
time, the combination of controllability and identifiability gramian is more ac- 
curate and has lower online times. The joint gramian supersedes in accuracy 
which is bought by a high offline time. 

The empirical cross gramian for MIMO systems completes the set empirical 
gramians for state reduction, while the joint gramian completes the body of 
parameter identification gramians. 
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Figure 5: Exemplary absolute error comparison of combined state and pa- 
rameter reduction for a linear system by balanced truncation (with empirical 
sensitivity gramian and empirical observability gramian), balanced truncation 
(with empirical identifiability gramian and empirical controllability gramian) 
and empirical joint gramian. 
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